GAUSSIAN BEAM METHODS FOR THE HELMHOLTZ EQUATION 

HAILIANG LIU\ JAMES RALSTON^ OLOF RUNBORG^ AND NICOLAY M. TANUSHEV* 

Abstract. In this work we construct Gaussian beam approximations to solutions of the high 
frequency Helmholtz equation with a locahzed source. Under the assumption of non-trapping 
fsf^ rays we show error estimates between the exact outgoing solution and Gaussian beams in terms 

,_^ of the wave number k, both for single beams and superposition of beams. The main result 

f^ is that the relative local L^ error in the beam approximations decay as k~^'^ independent of 

^Nj dimension and presence of caustics, for A'^-th order beams. 

^ 1. Introduction 

,__, In this article we are interested in the accuracy of Gaussian beam approximations to solutions 

-^ of the high frequency Helmholtz equation with a source term, 

^ (1) LnU =def Au + {iak + A;^)n^n = /, x G R'^. 

"^ Here k > is the wave number, assumed to be large, n{x) is the index of refraction and f{x; k) is 

^ a source function which in general also depends on k. We assume that both /(x; k) and n(x) — 1 

vanish for |x| > R. The nonnegative parameter a represents absorption. It is zero in the limit 
of zero absorption, where L^ solutions of (II]) become solutions satisfying the standard radiation 
condition. 

The Helmholtz equation (IT]) is widely used to model wave propagation problems in application 
areas like electromagnetics, geophysics and acoustics. Numerical simulation of Helmholtz be- 
^v^i' comes expensive when the frequency of the waves is high. In direct discretization methods a large 

^-H number of grid points is then needed to resolve the wave oscillations, and the computational cost 

■^ to maintain constant accuracy grows algebraically with the frequency. The Helmholtz equation 

^^ is typically even more difficult to handle in this regime than time-dependent wave equations, as 

" ' numerical discretizations lead to large indefinite and ill-conditioned linear systems of equations, 

• • for which it is difficult to find efficient preconditioners |12| . At sufficiently high frequencies direct 

simulations are not feasible. 

As an alternative one can use high frequency asymptotic models for wave propagation, such 
^ as geometrical optics |291 ITU ITT] , which is obtained when the frequency tends to infinity. The 

solution of the partial differential equation (PDE) is assumed to be of the form 

(2) u = ae''"'', 

where (p is the phase, and a is the amplitude of the solution. In the limit A; — )• oo the phase 
and amplitude are independent of the frequency and vary on a much coarser scale than the full 
wave solution. They can therefore be computed at a computational cost independent of the 
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frequency. However, a main drawback of geometrical optics is that the model breaks down at 
caustics, where rays concentrate and the predicted amplitude a becomes unbounded. 

Gaussian beams form another high frequency asymptotic model which is closely related to 
geometrical optics. However, unlike geometrical optics, the phase (j) is complex-valued, and 
there is no breakdown at caustics. The solution is still assumed to be of the form Q, but it 
is concentrated near a single ray of geometrical optics. To form such a solution, we first pick a 
ray and solve systems of ordinary differential equations along it to find the Taylor expansions of 
the phase and amplitude in variables transverse to the ray. Although the phase function is real- 
valued along the central ray, its imaginary part is chosen so that the solution decays exponentially 
away from the central ray, maintaining a Gaussian-shaped profile. For the simplest first order 
beams the phase <j) is a second order Taylor expansion, while the amplitude a is a zeroth order 
expansion. For wave equations one can use time as a parameter for the rays, and the expressions 
for the phase and amplitude are 

(3) 0(i, y) = Mt) + {y- x{t)) ■ p{t) + 2 (y - <t)) ■ M{t){y - x{t)), a{t, y) = ao{t) 

where x{t) is the geometrical optics ray, p{t) is the direction of the ray and the second derivative 
matrix M{t) encodes the width and curvature of the beam; M has a positive definite imaginary 
part which ensures the beam has a Gaussian shape. In the Helmholtz case, since there is no 
longer a distinguished variable with level sets transverse to the rays, one uses Taylor expansion 
in the plane orthogonal to the ray direction. Higher order beams are constructed through higher 
order Taylor expansions in ([3]). 

The existence of Gaussian beam solutions to the wave equation has been known since sometime 
in the 1960's, first in connection with lasers, see Babic and Buldyrev [2]. Later, they were used 
in the analysis of propagation of singularities in PDEs by Hormander [21] and Ralston [39] . In 
the context of the Schrodinger equation first order beams correspond to classical coherent states. 
Higher order versions of these have been introduced to approximate the Schrodinger equation 
in quantum chemistry by e.g. Heller [16], Hagedorn [H], Herman and Kluk |17j . 

More general high frequency solutions that are not necessarily concentrated on a single ray 
can be described by superpositions of Gaussian beams. This idea was first introduced by Babic 
and Pankratova in [3j and was later proposed as a method for approximating wave propagation 
by Popov in [41j. Letting the beam parameters depend on their initial location z, such that 
X = x{t]z), p = p{t\z) etc., and a = a{t,y;z), (p = (p{t,y;z), the approximate solution for an 
initial value problem can be expressed with the superposition integral 

(4) n(t, y) = ( A j ' j^ ait, y; z)e*'=*(*'^^^)dz , 

where Kq is a compact subset of M'^. 

It should be mentioned that there are other related Gaussian beam like approximations. In 
the thawed Gaussian approximation [15] the phase cj) is always a second order polynomial. Higher 
order is obtained by instead taking a higher order polynomial in the amplitude, to correct also for 
errors in the phase. Frozen Gaussian approximations [161 [T7] also use a second order polynomial 
for the phase (/>, but with a fixed size of the second derivative (M(t)=constant). Single frozen 
Gaussians are therefore not asymptotic solutions to the wave equation. However, superpositions 
of frozen Gaussians are and they can be thought of as an efficient linear basis for the wave 
equation. 

Numerical methods based on Gaussian beam superpositions go back to the 1980's with work 
by Popov, [4H I27j . Cerveny [10| and Klimes |30j for high frequency waves and e.g. Heller, 
Herman, Kluk |16l [T7| in quantum chemistry. In the past decade there was a renewed interest in 



such methods for waves following their successful use in seismic imaging and oil exploration by 
Hill |181ll9j. Development of new beam based methods are now the subject of intense interest in 
the numerical analysis community and the methods are being applied in a host of applications, 
from the original geophysical applications to gravity waves [l6], the semiclassical Schrodinger 
equation [13^ [2^ [3T] , and acoustic waves |45j . See also the survey of Gaussian beam methods 
in |23j . Individual beams are normally computed in a Lagrangian fashion by solving ODEs 
along the central rays. The superposition is then replaced by a discrete summation of beams. 
There are also more recent numerical techniques based on Eulerian formulations of the problem 
[62\ [2^ [221 EU I32]- In these methods a PDE is derived for the parameters in the beams, i.e. the 
quantities in the ODEs. This is coupled with a level-set PDE for the ray dynamics. With the 
Eulerian formulation the result is no longer a superposition of asymptotic solutions to the wave 
equation! For superpositions over subdomains moving with the Hamiltonian flow, it was shown 
directly in [341 [35j that they are asymptotic solutions without reference to standard Gaussian 
beams. Numerical approaches for treating general high frequency initial data for superposition 
over physical space were considered in [Ul [1] for the wave equation. 

In this paper we study the accuracy in terms of k of Gaussian beams and superpositions of 
Gaussian beams for the Helmholtz equation (II| . This would give a rigorous foundation for beam 
based numerical methods used to solve the Helmholtz equation in the high frequency regime. 
In the time-dependent case several such error estimates have been derived in recent years: for 
the initial data |45j , for scalar hyperbolic equations and the Schrodinger equation [Ml |35l [36] , 
for frozen Gaussians |43| [33] and for the acoustic wave equation with superpositions in phase 
space [7] . The general result is that the error between the exact solution and the Gaussian beam 
approximation decays as k~ ''^ for A'^-th order beams in the appropriate Sobolev norm. There 
are, however, no rigorous error estimates of this type available for the Helmholtz equation. What 
is known is how well the beams asymptotically satisfy the equation, i.e. the size of L„u for a 
single beam. Let us also mention an estimate of the Taylor expansion error away from caustics, 
[38]. 

The analysis of Gaussian beam superpositions for Helmholtz presents a few new challenges 
compared to the time-dependent case. First, it must be clarified precisely how beams are 
generated by the source function and how the Gaussian beam approximation is extended to 
infinity. This is done in §2 and §3 for a compactly supported source function that concentrates 
on a co-dimension one manifold. Second, additional assumptions on the index of refraction n{x) 
are needed to get a well-posed problem with ^-independent solution estimates and a well-behaved 
Gaussian beam approximation at infinity. The conditions we use are that n{x) is non-trapping 
and that there is an R for which n{x) is constant when \x\ > R. 

In §4 we consider the difference between the Gaussian beam approximation and the exact 
solution to the radiation problem with the corresponding source function. Here we are interested 
in behavior of the local L^ norm \\ugb — '^\\l^{\x\<r) as A: — >• oo. This depends on the well- 
posedness of the radiation problem. There are a variety of estimates that apply here [iQl [8], but 
the Laplace-transform based estimates of Vainberg |48l H9] suffice for our purposes. In §5 we 
compare the Gaussian beam approximation with the result of stationary phase expansion of the 
exact solution in a simple example. 

Sections §6 and §7 are devoted to superpositions of beams with fundamental source terms. 



Our main result is Theorem 6.1 where we are able to show that the error between superposition 
of A'^-th order beams and the exact outgoing solution decays as k~ ''^ independent of dimension 
and presence of caustics. This is consistent with the optimal results of |36j in the time-dependent 
setting. Finally, §7 gives an example of how beams can be constructed for more general source 
functions. 



4 HAILIANG LIU\ JAMES RALSTON^, OLOF RUNBORG^*, AND NICOLAY M. TANUSHEV 

2. Construction of Gaussian beams 

In this section we construct the Gaussian beam solutions for M when / is compactly sup- 
ported on a co-dimension one manifold. This construction has become standard (see, for exam- 
ple, [39j or [28j ) and we review some details here which will be used later. The form of the beam 
solutions is 

(5) u{x; k) = e*'='^(^)(ao(x) + ai{x)k-^ + ■■■ + ae{x)k-^). 

Each beam concentrates on a geometrical optics ray 7 = {x{s) : s £ M}, which is the spatial part 
of the bicharacteristics {x[s),p{s)) defined by the flow for the Hamiltonian H{x,p) = |pp — n^(x) 

(6) x = 2p, ■p = —Vxn^{x). 

We assume that there is a number R > such that the (smooth) index of refraction satisfies 
n{x) = 1 when |x| > R and that the source function / is compactly supported in {|a;| < R}. Here 
we also restrict the construction of the Gaussian beam solution to the larger region \x\ < 6R. 
The essential additional hypothesis for our construction is that the index of refraction does 
not lead to trapped rays. The precise non-trapping condition is that there is an L such that 
\x{L)\ > 2R for all solutions with |x(0)| < R and H{x{0),p{0)) = 0. Note that this implies that 
\x{s)\ > 2R for s > L since rays are straight lines when n{x) = 1. 
Applying L^ in ([I]) to ([s]) we have 

e 

(7) L„^ = e^'=<^ J;c,■(x)fc-^ 

where 

c_2 = (n^ - |V^(/)p)ao =def E{x)ao, 

c_i = ian^ao + V^,. • (aoVa,.0) + Vxao ■ ^ x(p + Eai, 

Cj = ian^aj+i + Vx ■ {aj+i\7x4>) + ^xaj+i ■ ^x4> + Euj+i + AxQj, j = 0,1,. ..,£. 

ODEs for S{s) = (j){x{s)), M{s) = D'^(j){x{s)) and Aq{s) = ao(x(s)) arise from requiring that 
c_2 vanishes to third order on the ray x{s), and that c_i vanishes to first order on the ray. It 
leads to the equations 

(8) S = 2n^{x{s)) M = D^{n'^){x{s))-2M^ Ao = -tT{M{s))Ao - an^{x{s))Ao. 

This amounts to constructing a "first order" beam. Higher order beams can be constructed by 
requiring c_2 vanishes to higher order on 7. Then one can require that the Cj's with j > —2 also 
vanish to higher order, and obtain a recursive set of linear equations for the partial derivatives 
oi ao,ai, . . . ,ai. More precisely, for an A^-th order beam £ = [A^/2] — 1 in (Isj) and Cj{x) should 
vanish to order N — 2j — 2 when —2 < j < £ — 1- 

For initial data, we let 5(0) = and choose M(0) so that 

(9) M(0) = M(O)'^, M(0)x(0) = p(0), Im{Af (0)} is positive definite on x(0)-^. 

Then for all s the matrix M{s) inherits the properties of M(0): M{s)x{s) = p{s), M{s) = 
M{s)^ , and Im{M(s)} is positive definite on the orthogonal complement of x{s), see |39j . For 
the amplitude we take ^o(O) = 1. We can solve the ODE for Aq explicitly, and obtain 



Aq{s) = exp f - / {an^{x{T)) + tr(Af ))(ir 

The phase </) in ([s]) can be any function satisfying (j){x{s)) = S{s), V(f){x{s)) = p{s) and 
D'^(f){x{s)) = M{s). However, to write down such a function we need to have s as a function of 




p>0 



Figure 1. Notation for the source in two dimensions. The gray area indicates ^{i])- 



X. Since we have x{s) ^ 0, x(s) traces a smooth curve 7 in W^, and the non-trapping hypothesis 
imphes that this curve is a straight hne when |s| > L. We let 

r2(r/) = {x : \x\ < 6R and |x — 7I < r]}, 

be the tubular neighborhood of 7 with radius r] in the ball {|x| < 6R}. By choosing t] small 
enough, we can uniquely define s = s{x) for all x G ^{r]) such that x{s) is the closest point on 
7 to X, provided 7 has no self-intersections. We then define the phase function and amplitude 
A on Q for first order beams by 

(10) 0(x) = S{s) +p{s) ■ {x - x{s)) + -{x - x{s)) ■ M{s){x - x{s)), A{x) = Ao{s), 

with s = s{x). Note that s{x) is constant on planes orthogonal to 7 intersected with il(7?). 
Since 7 can have only finitely many self-intersections, we can cut 7 into segments without self- 
intersections, and define s{x) on a tubular neighborhood each segment, ignoring the endpoints. 
For this reason self-intersections will not create difficulties, and without loss of generality we 
will assume that 7 has no self-intersections in what follows. The construction of the Gaussian 
beam phase and amplitude for higher order beams is carried out in a similar way |39j. 

2.1. Source. To introduce the source functions that we will consider in this article let p be a 
function such that \Vp\ = 1 on {x : p{x) = 0}, and define S to be the hypersurface {x : p{x) = 0}. 
Given xq S H, we let {x{s),p{s)) be the solution of (6) with (2;(0),p(0)) = (xo,n(xo)Vyo(xo)). 
Since we assume no trapped rays and n{x) = 1 when |x| > R, x{s) and p{s) are defined for 
s E M, and we set 7 = {x(s), s E M}. Then we can assume that s{x) is defined on the tubular 
neighborhood 0(77) of 7 as above (assuming no self- intersections) . We begin with a beam u(x, k) 



concentrated on 7, and defined on ri(r/). If u is first order, we can define it by (10). Then we 
define u^ to be the restriction of u to {x : p{x) > 0}. In order to have a source term which is 
a multiple of 5{p), we need a second beam u~{x, k) defined on {x : p{x) < 0} which is equal to 
u~^ on S for all k. Hence, writing u^{x,k) = A~^{x,k)e^'"l' ^^^ and u~{x,k) = A~{x,k)e^^'^ ^^\ 
we must have 0^ = (f)~ and A'^ = A~ on S. Those requirements and Cj = 0, j = —2,...,£ 
at xq determine the Taylor series in the transverse variables at xq for (/)^ and A^ . To see this 
suppose that u~ is going to be a beam of order A^ and that the coordinates on ^{rj) are given 
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by {s, y) where s = s{x) and y = (yi, . . . , yd-i) is transversal. Then, provided rj is chosen small 
enough, S is given by s = cr(y) with o"(0) = and V(t(0) = 0. To determine the Taylor series 
in y for (j)~{s,y) at s = one differentiates the equation (j)~ {a{y),y) = 4)'^{a{y),y) with respect 
to y and evaluates at y = 0. When partial derivatives of (j)~ with respect to s appear in this 
calculation, they are determined by the requirement that c_2 vanishes on x{s) to order A^ + 2. 
The Taylor series for A~ in the transverse variables at xq is determined in the same way from 
A^{a{y), y, k) = A^{a{y), y, k) for all k. To construct u~ , we use those Taylor series as data at 
s = in solving the equations Cj = 0, j = —2, . . . ,i along x{s). Since for an A^-th order beam 
we only require that Cj vanishes on x{s) to order N — 2j — 2, we can still require that (/>"'" = (/>~ 
and A~^ = A~ exactly at points on S. Extending n"*" to be zero in {x : p{x) < 0} and u~ to be 
zero in {x : p[x) > 0}, we define ugb = u^ + u~ . Then we have, setting A = A'^ = A" on S, 



(11) LnUGB 



ik 



di 



di 



A + 



dA+ 
dv 



dA- 
dv 



e'"* S{p) + fGB=defgoS{p) + fGB, 



where v{x) = 'Vp{x), the unit normal to S. We consider the singular part of LnUGB in (H)) i-6- 
go5{p), to be the source term and Jgb to be the error from the Gaussian beam construction. 
Note that 



(12) 



fcB 



Jk(P+{x) V^ 



i=-2 



C. [X 



\k ^ + 



Jk(p (x) 



Yl S' (^)^~ 



where the c^{x) are extended to be zero when p{x) < and the c- (x) are extended to be 
zero when p{x) > 0. For first order beams i = and <\8h implies c_2{x) and c_^(x) are 
0{\x — x(s(x))p) and 0{\x — x(s(x))|) respectively. Finally we restrict the support of ugb to 
^iv) by multiplying it by a smooth cutoff function supported in r2(ry) which is identically one on 
the smaller neighborhood 0(77/2). The cutoff function modifies A^, and Jgb, outside Q,{rj/2), 
but its contribution to (11) is exponentially small in k (see |36)). and we will disregard it from 
here on. 



2.2. Estimate of Jgb- From the non-trapping condition, it follows that the length of a ray 
inside ^{rj) is bounded independently of starting point in |x| < R. By construction, c^ (x) is 
bounded and 



(13) 



^?(^)= E d%^{x){x-x{s)f, 

m=N-2j-2 



-2,...,i-l, 



where dg (x) are bounded on il(??). Hence, 



"/3,i^ 



\cf{x)\ <Cj\x- x(s)|^-2j'-2, X E 0(77). 



Choosing t] sufficiently small, the construction also ensures that 

(14) Im{(/)=^}(x) > c|x -x(s)p, x£Q{r]), 
see [36]. From the bound 

(15) sPe""^' < C7pa-P/2g-asV2^ ^^ ^ {p/ef/'^, 




R 2B. 3R m 5R GR 



Figure 2. The cut off functions ?/3(x) and rj^ix) 



with p = N — 2j — 2, a = kc and s = \x — x{s)\ we then get for x £ ^2(77), 

e e 

|/GB(a;)| <e-'=^"^^<^*>(") Y^ |c±(x)|A;-^' < e-^^l"-"(^)l' Yl Cj\x - x{s)\^-^^-^k-^ 

i=-2 j=~2 

e 

(16) < (J^-kc\x-x{s)\y2 y^ j^-N/2+j+lj^-j ^ ^^-kc\x-x(s)\y2j^-N/2+l_ 

J=-2 

We note that the constant is uniform in |x| < 6R and in particular for first order beams fcB 
willbeO(A:i/2e-fc£|^-^WP). 

3. Extension of Gaussian beam solutions to infinity 

In this section we extend ugb{x) defined on |3;| < QR to an outgoing solution ugb{x) in M . 
For estimates on the validity of the approximation it is essential to do this so that 

IgB =def LnUGB - go^ip), 

is supported in |a;| < QR and is o{k). 

The main step in the extension is a simplified version of the procedure used in |37] . Let G\{x) 
be the Green's function for the Helmholtz operator A + A^, where A may be complex valued. 
When a > 0, define 



(17) 



ka '■= yk'^ + ika. 



Then Li = A + iak + /c^ = A + /c^, and Gk^ is uniquely determined when q > as the inverse 
of the self-adjoint operator Li; for q = it can be defined either as limQ,4,o G^fca or by radiation 
conditions. In the case d = 3, 



Gfcjx) = -(4^)-i 



s^/'-al^l 



\X\ 



1 \x\<{a-l)R 
\x\ > aR 



To extend uqb we introduce the cutoff function r]a{x) in C°°(M ) with parameter a > 1: 

r]a{x) 
(see Figure [2]) and define 

(18) UGB = mix)uGBix) + I Gk^{x - y)r?5(y)L„[(l - m{y))uGB{y)]dy. 
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We also assume that R is chosen large enough such that the support of goS{p) C S n ^{rj) is 
inside {|x| < R}. 

Consider first LnUcB in the region {|x| > R}. Since Ln = Li as well as go5{p) = in this 
region and r/5 = 1 on the support of r/3, 

fcBix) = LnUGsix) = 7]5{x)Ln['n3{x)uGBix)] + r/5(x)L„[(l - T]3{x))ugb{x)] 
= 'q^{x)LnUGB= V5ix)fGB{x). 

Since % is supported on |2;| < 5R, it follows that /gb vanishes for \x\ > 5R. 



Consider next the region{|rE| < R} and let v = ugb — Vs'^^GB, i-e. the integral term in (18). 
Since, 7/3 = 1 on |x| < R, we have in this region 

UGB - UGB = V, fcB - fcB = LnV. 

In view of the estimate of Jgb it now suffices to show that for |x| < R, dxV decays rapidly when 
/c —7- 00, for all multi-indices, |/3| < 2. 

By the definition of the two cut-off functions, we have for |x| < i? 

v{x) = / Gk^{x - y)r]5(y)Ln[il - m{y))uGB{y)]dy 

Gk^x - y)v5{y)Li[{i - m(.y))uGB{y)]dy. 



'2R<\y\<5R 
The fundamental solution Gu has the form 

d-3 d~3 

where w and its derivatives in x are bounded by |/ca| 2 < Ck 2 on compact subsets of |a;| > it!, 
see Appendix. Since n{x) = 1 for |x| > R, in that region x{s) is a straight line and Vx4'^{x{s)) 
is a constant unit vector. Since x{s) is going out of |a;| < R when it crosses \x\ = R, at x{s) = y 
with \y\ > 2R the phases in ugb satisfy V^i?i>^(x(s)) • y > cos(7r/6)|y|. Likewise when \x\ < R 
and \y\ > 2R, {y — x) ■ y > \y\\y — x\ cos(7r/6) (see Figure [s]). The form of ugb (see ^) gives 
the integrand in (18) the form e^^^b{y, k) with ip{y) = 4'^{y) + {ka/k)\x — y\ and b smooth in y, 

d-3 

bounded together with its derivatives by Ck 2 . Note that 

k \y — x\ 

The preceding remarks show that, when |x| < R and k large, Vytp does not vanish on the 
support of the integrand in (18). Hence we can use the identity 

and integrate by parts to show that v and its derivatives are order k~"^ for any m. 
This completes the verification of the extension. We have shown that 

(19) fGB{x)=m{x)[fGB(.x)+r{x)], \\r\\Lmx\<m)=0{k~'^). 

Hence, the size of /gb is of the same order as the size of /gB; which is 0{k^ ''^^^e~ '^~^^-^'' ). 
Moreover, 

(20) \\uGB - ugbUl^IxKR) = Oik'""), 

for any m. Note that, since ugb is represented by Ga,k for |x| large, it is square-integrable 
(a > 0) or outgoing (a = 0). 




Figure 3. Maximum angle. 



4. The Error Estimate for ugb 

In this section we will use an estimate showing that the radiation problem is well-posed due 
to Vainberg [l8] and [l9] . This will give estimates on the accuracy of ugb as an approximation 
to the exact solution in the region |x| < R. Vainberg starts with the initial value problem for 
wave equation in M^ x Mj 



vtt 



n 



'-Av = 0, v{0) = 0, vt{0) 



-n ^5 



and takes the Fourier-Laplace transform 
(21) u{x,k) 

to get the solution to 



jxt 



v{t,x)dt 



Am + A^n^n 



satisfying radiation conditions. Taking advantage of finite propagation speed, and the prop- 
agation of singularities to infinity, he can estimate u on bounded regions from the integral 



representation (21), when g has bounded support and the nontrapping condition holds. In the 
notation of [48j, u = [TZx]{n~'^g), where Tl\ is the operator 

This is defined for complex A as the analytic continuation of T^a restricted to the space HJ^ with 
range in the space -ff™(|x| < b). The estimates take the following form: there are constants C 
and T such that 



(22) 



II^A5lU+2-,,(6)<C|A|i-^-e^|i-^l| 



\9\ 



„a, < J < 3. 



Here the norms are standard Sobolev norms on H^{M. ), the closure of C^(|x| < a) in || • ||m, 
and H"^{\x\ < b). One can assume that b < a. The admissible set of A here is the set 



U, 



Cl,C2 



{A G C : |Im A| < d log |Re A| - C2} 



for some ci, C2 > 0. If d is even, then one has to add the condition 

-7r/2 <argA< 37r/2. 
This is Theorem 3 for d odd and Theorem 4 for d even in 
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Here we will apply (22) with g = n "^fcB, cl = 6R, b = R and A= /cq, G C with ka defined in 
(17). This makes v?TZka9 = 'Ugb — ue, where ue is the exact solution to the radiation problem 



(TFwith / = goS{p) defined in (11). Taking m = and j = 2, we have 
(23) \\uGB-UE\\m\.\<R)<C\k^\-'e^\'^''-\\\fc 

Note that \ka\ = k {l + {a/kfY^^ and 

k ( 



Im kn 



M^^+i-m'" 



1 



1/2 



Re kr, 



V2\ 



1 



[cik?r-^i'"' 



Hence |Im ka\ < C , ka & f^ci,c2 for some ci, C2 > and \ka\ > k, so 
(24) II--- --11 ^ ^ r^iu-h ~ 



\UGB - Ue\\l-^{\x\<R) < C\k\ ||/gb||l2 + \\UGB - ""Gsl |l2(|x|<R) > 



uniformly in terms of a. The estimates in ( |19| ) and (20) ensure that 

\\UGB - Ue\\l'2{\x\<R) <C\k\~ ||/gb||l2(|^|<5^). 



(25) 



We observe here that since (19) and (20) hold uniformly for all beam starting points xq G S 



the estimate (25) will also hold for linear superpositions of beams, which we will discuss further 
below, see (|35[). Moreover, from (16) and the estimate (45) derived below, we obtain 



ll/GB|li.(|,.|<5iJ)<CA;-^+' 



g-2fcc|x-x(s)|2^^ ^ ^^-Ar+2+(l-d)/2 



Q.{-n) 



This finally shows that for a single beam ugb-, 



d-l 



o-d 



Note that the factor k '^<' corresponds to the size of the L^ norm of the beam itself in d di- 
mensions, ||wgb||^2(|^|^^) ~ k^"''^, showing that the relative error of the beam is bounded by 

5. An Example 
Using the notation x = {xi,x') = (xi, X2, a^s), the outgoing solution to 

An + A;\ = 2iA;e-'=l^'l'/25(xi) 
is given by 

^ik\x-{0,y')\-k\y'\y2 



(26) 



m(x, k) 



-2ik 



-dy. 



4vr J^2 |a;-(0,y')l 

In this section we compare the approximation that one gets by using the method of stationary 
phase on this integral to the approximation given by ugb- The stationary phase approximation 
is not uniform in x' , and for x' 7^ it simply gives u[xi,x',k) = 0{k~ ) for all A^. However, 
when x' = 0, it gives ugb{xi^{^). 

The procedure for constructing u^ given earlier with the source 2zfee~'^ ' "(5(xi), gives x{s) = 
(2s, 0,0), p{s) = (1,0,0), S{s) = 2s, M{s) = j^P and A{s) = {I + 2is)-\ where P is the 
orthogonal projection on ej-. For u~ one gets the same results with s replaced by —s and p{s) 
replaced by —p{s). The definition of s(x) gives s{x) = |xi|/2, and we have 



(27) 



UGB{x,k) = {l + i\xi\) ^e*'"^, where 



\xi\ + 



/|2 



2(l + i|xi| 
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To apply stationary phase to (26) assume that xi ^ 0. Then the phase is given by 'ijj{x,y') 
\x-{0,y')\+i\y'\'^/2 and 



y'-x' 



^ \x-{0,y')\ 
That vanishes and is real only when y' = x' = 0. Then one has 



1py'y'\x'=y'=0 

The stationary phase lemma (P^) gives 

(28) 

Since 



1 

kil 



+ i]I' 



'2x2- 



n(xi,0) = (det(-i^,v(xi)))-V2 ( -^^-^ + 0(1) ) . 



k 



A'K \x\ 



<let{-iiiyiyi{xi)) = [ — — + 1 ) , 

K'l 



and the choice of square root leads to 



+ 1 



\xi\ 



-1/2 



+ 1 



Fl 



-1 



one sees that the leading term in (28) is exactly (27) 



6. Error Estimates for Superpositions 



Given a point z G S, we relabel the primitive source term qq in (11) as 
(29) g{x, z, k) = [ikUx) + C2(x)]e-'=l--^l'/2<^(p), 

where Qj G C^ and Ci(^) = 1 on a neighborhood oi x = z. Denoting the resulting beam as 



ugb{x', z), the error estimate (24) is uniform in z as long as z remains in a compact subset of 



|x| < R, for instance \z\ < R/2. If we let z range over E, we can form 

(30) g{x,k)dip) 

and 

(31) 



^ N (d-l)/2 . 

— j / gix,z,k)h{z)dA^ 



u{x) 



211 J 



{d-l)/2 



ugb{x; z)h{z)dAz 



is an approximation to the exact solution for the source g{x, k)5{p) satisfying the estimate (24). 



We now state the main result of error estimates for superposition (31). 



Theorem 6.1. Assume thatn{x) is smooth, non-trapping, positive and equal to 1 when \x\ > R. 



Let ue he the exact solution to (Tn) with the source f = g{x, k)5{p) in (30), and u he the Gaussian 



heam superposition defined in \31y hased on N-th order heams. We then have the following 
estimate 

(32) \\u-UE\\L^i^\.\<R)<Ck~'"\ 

where C is independent of k hut may depend on R. 
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In order to simplify the notation, we specify p{x) = xi and y = (0, z) for z G S C M ^. The 
superposition thus can be written as 

^ X (d-l)/2 r. 

— j / UGB{x;z)h{z)dz, 



u(x) 



(33) 

and the residual 

(34) LnU - LnUE = f {x) 



— j / fGB{x;z)h{z)dz. 



By the definition of ue and the source g{x, k)6{p), the residual / contains only regular terms. 
We can therefore extend the superposition u to u in the same way as in §3, and define / = 
LnU — LnUE- As observed above, (19) and (20) hold uniformly for all z G S, and the same steps 
as in §4 therefore lead to an estimate corresponding to (25), namely 

(35) ||n - ue\\l^(\x\<r) < Cfc~"^ll/llL2(|^.|<5ij). 

We let x{s; z) be the ray originating in z, x(0, z) = z and we denote by ^{rj; z) the corresponding 
tubular neighborhood of radius r/, in the ball {|a;| < 5R}. By choosing 77 > sufficiently small, 
we can thus ensure that s = s{x; z) is well defined on Q{r];z). In what follows we denote 
x{s{x, z); z) by 7 or 7(x; z). Moreover, we introduce the cutoff function ^.^(x) G C^iW^) as 

1 for |x| < 77/2, 
for |x| > r/. 



(36) 



Qr)[ 



> and Qn{x) = \ p. 



such that gj^{x — ^{x; z)) is supported on $7(r/; z) and is identically one on 0(r?/2; z). The form 
(12) of fcBix; z) will then be 



fGB{x;z) 



^ik(p^ (x;z) 



E 



cj{x;z)k ^ +e 



ik(j> {x;z) 



Y, c-{x; z)k-^ Br, (x - 7) + 0(fc-°°) 



with bounds 



Y, k^''e'^^-^''''Uaix] z)ix - -ff^Br^ (x - 7) + 0(A:-°°), 



|/3„|<iV + 2, 2j„<2-iV + |/?„|. 



The sum over a is finite, da involves the functions do ■ in (^13| and (pa is either cf)^ oi cj) . 
Moreover, 0{k~°^) indicates terms exponentially small in \/k. After neglecting these terms and 
using ( 34 ) it follows that we can bound the L^ norm of / by 

2 

-JV+I/Sq ., , r. 

2 e"^'^"' da{x - -iY°' Br^hdz 



L2(|x|<5fi) 






L2(|x|<5iJ) 



Ia[x^ z, z')dzdz'dx, 



ck"-''' 

^ J\x\<5RJT: JY, 
where the terms la are of the form 

Iaix,z,z') = k'^+\^\e'^^^'^'''''^ g{x; z')'g{x^ 

X (x- ^f [x - 7) Bn {x - 7) £»^ [x - 7 ) , 1^1 < 3. 
Here ^(x; z) = da(x; z)h{z) and 
(37) il){x,z,z') := (t){x;z')- (t){x;z), 
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with (]) being either of cj) . The function g and its derivatives are bounded, 



(38) 



sup \d^g{x;z)\ < Cx, 

ze'S,xen{'ri;z) 



for any |A| > 0. 

Let Xji^j ^j z') ^ C°° be a partition of unity such that 



Xi{x,z,z') 



1, when |7(x, z) — 7(0::, z')| > 6*12; — z'l, 



0, when 17(2;, z) — j{x, z')\ < \0\z — z'\ 
and Xi + X2 = 1- Moreover, let 

h = Xl{x, Z, z')Ia{x, Z, z'), h = X2{x, Z, z')Ia{x, Z, z'), 

SO that Ia{x, z, z') = Ii + I2. 

The rest of this section is dedicated to estabUshing the following inequality 



(39) 



L 



\x\<^RJY, JT, 



Ij{x, z, z')dxdzdz' 



<Ck 



2-d 



for j = 1, 2. With this estimate we have ||/||l2(|^|<5^) < Ck^^'^''^, which together with (35) lead 
to the desired estimate ( 32 ) . 

A key ingredient in establishing estimate ( 39 ) is a slight generalization of the non-squeezing 
lemma obtained in [36] . It says that the distance in phase space between two smooth Hamiltonian 
trajectories at two parameter values s that depends smoothly on the initial position z, will not 
shrink from its initial distance, even in the presence of caustics. The lemma is as follows: 

Lemma 6.2 (Non-squeezing lemma) . Let X = (x{s; z),p{s;z)) be the bi- characteristics starting 
from z E T, with T, bounded. Assume that p{0; z) £ C^(S) is perpendicular to E for all z, that 
\p(0;z)\ = n{z) and that inf ^ n(2;) = no > 0. Let S{z) be a Lipschitz continuous function on S 
with Lipschitz constant Sq. Then, there exist positive constants ci and C2 depending on L, Sq 
and no, such that 

(40) ci\z - z'\ < \p{S{z); z) - p{S{z'); z')\ + \x{S{z); z) - x{S{z'); z')\ < C2\z - z'\ , 

forallz,z'G^ and\S{z)\,\S{z')\ < L. 



Proof. With the assumptions given here, the non-squeezing lemma proved in |36| states that 
there are positive constants < di < d2 such that 



(41) 



di\z — z'\ < \p{s; z) — p{s; z')\ + |x(s; z) — x{s; z')\ < d2\z — z'\ 



for all z,z' €z T, and |s| < L, i.e. essentially the case S{z) = constant. Since the Hamiltonian 
for the flow (p| is regular for all p, x, and the initial data p{0; z) is C^(S), the derivatives <9°^x 
and d"^p with \a\ < 2 are all bounded on [—L,L] x S by a constant M. Then, for the right 
inequality in (|40|), we have 



\p{S{z);z) - p{S{z');z')\ + \x{S{z);z) - x{S{z'); z')\ 

< \p{S{z);z)-p{S{z'y,z)\ + \p{S{z'y,z) -p(S{zy,z')\ 

+ \x{Siz);z) - x{S{z'y,z)\ + \x{S{z'y,z) - x(5(z'); z')\ 

< 2M\S{z) - S{z'y + d2\z - z'\ < {2MSo + d2)\z - z'\ =: C2\z 
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by ( |41[ ) and the Lipschitz continuity of S{z). For the left inequahty in (40), 
(42) IxiSiz); z) - x{S{z')-z')\ + \p{S{z)-z) - p{S{z'); z')\ 

> \p{S{z)-z) -p{S{z)-z')\ - \p{S{z)-z')-p{S{z')-z')\ 
+ \x(S{z);z) - x{S{z);z')\ - \x{S{z); z) - x{S{z'); z')\ 

>di\z- z'\ - \p{S{z);z') - p{S{z');z)\ - \x{S{z); z) - x[S{z)\z')\ 

>di\z-^\ -2M|5(z)-5(z')|, 



where we again used (41). Next we will estimate \S{z) — S{z'y\ using \x{S{z')\z') — x{S{z)\z)\. 
From Taylor expansion of x around z, and the fact that Xs = 2p, we have 

x{S{z')-z') - x{S{z)-z) = 2p{S{z)- z){S{z) - S{z')) + D,x{S{z)-z){z' - z) + R{z, z'), 

where 

(43) \R{z, z')\<M {\S{z) - S{z')\^ + \z- z'\^) < M(l + Si 



2m„ /|2 



0^ 



z — z 



Moreover, 

-p{s; zY'Dzx{s] z) = ps{s; z)'^D;,x{s; z) + p{s; z^D^Xsis] z) 



ds 

= -Vxn^{x{s; z))'^Dzx{s] z) + 2p{s] z)^D-^p{s] z) 

= V,H{x{s; z),p{s; z)) = V,H{x{0; z),p{0; z)) = 0, 

by the choice of data at s = 0. Therefore, since p{0; z) is orthogonal to S and Xz{0;z) are 
tangent vectors to S, we have p{s; z) Dzx{s; z) = for all s and 

(44) \x{S{z);z) - x{S{z'y, z')\ > 2\p{S{zy, z)\\S{z) - S{z')\ - \R\ > 2no\S{z) - S{z')\ - \R\. 



Together (|42j), (|43|) and (|44|) now give 

\x{S{zy, z) - x(S{z'y z')\ + \p{S[z)- z) - p{S[z'); z')\ 

> d,\z - z'\ - —\x{S{zyz) - x{S{z'y,z')\ - ^^^"ilt^i^ _ z'W 
no no 

which implies 

\x{S{z); z) — x{S{z')\ z'y + \p{S{z); z) — p{S{z'); z')| > di\z — z'\ (l — m\z — z'|) , 

with m = M^(l + S'o)/(no(ii) and di = di/{l + M/rio). The lemma is thus proved for \z — z'\ < 
l/2m with ci = di/2. On the other hand, if jz — z'l > l/2m there is a number c{m) such that 

inf |p(s; z) — p{s'] z')\ + \x{s; z) — x{s] z'y =: c{m) > 0, 

2,^'eS, \z-z'\>l/2m 
\s\<L, \s'\<L 

by the uniqueness of solutions to the Hamiltonian system. Hence, in particular, for \z — z'\ > 
l/2m, 

\xisizyz) - xisiz'yz'y + Hsizyz) - p{s{z'y z'y > c{m) > ^|z - /|, 

where A = sup^ ^i^y. I-^ — -z'I < oo is the diameter of the bounded set S. This proves the lemma 
with ci = min((ii/2, c(m)/A). D 



We now prepare some main estimates for proving (39). 
Lemma 6.3 (Phase estimates). Let t] be small and x G D(ri,z,z') where 

D{r], z, z') = n{ri, z) D n{ri, z'). 
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• For all z,z' ^Ti and sufficiently small r], there exists a constant 5 independent of k such 
that 

QiIj [x , z , z') > 6 |x — 7I +I2; — 7'! 

• For \^{x; z) — ^{x; z')\ < 9\z — z'\, 

\V,^l^ix,z,z')\>C{e,r])\z-z'\ , 
where C{6, rj) is independent of x and positive if 9 and r] are sufficiently small. 



Proof. The first result follows directly from (14). For the second result, we proceed to obtain 



\Vx'4'{x,z,z')\ > \^Vxip{x,z,z')\ 

= l^V^cPix; z') - 3ftV,(/.(x; z)|, [h := SRV,^} 

= h{-f';z')-hir,z)+h{r,z')-h{^';z') 

+ h{x; z) - h{'-f; z) + /i(7, z) - h{x, z) . 

For the function z 1— )• s{x; z) we can find a Lipschitz constant that is uniform in x. Recalling 



that 7 = x{s{x; z); z) and 7' = x(s(x; z'); z') we can therefore use (40) in Lemma 6.2 for the first 
pair, and obtain 

1/1(7; z') -h{r,z)\ = \p{s{x;z)]z) -p{s{x;z');z)\ > ci\z - z'\ - |7-7'l- 

The second pair |/i(7; z') — h{j'; z')\ is bounded by C1I7— 7']. Then, by the Fundamental Theorem 
of Calculus, for x £ D(rj, z, z'), the remaining terms are 



/ [D^(I){tx + (1 - r)7; z') - D'^cpirx + (1 - T)r, z)] (x - 7)dr 
Jo 



< C\z — z'\\x — 7I < C2'i]\z — z'\ 



Using these estimates for the case I7 — 7'] <0\z — z'\ we then obtain 

|Vx^(x, z,z)\'>c\\z- z\ - I7 - 7'! - C1I7 - 7'! - C2r]\z - z'\ 
>ci\z-z\-{l + Ci)e\z - z\ - C2v\z - z'\ 
=:C{e,r])\z-z'\ , 

where C{9, rj) is positive if 9 and r] are small enough. 



D 



6.1. Estimate of /i. We start by looking at /i which corresponds to the non-caustic region of 
the solution. We have 



/i(x, z, z)dzdzdx 

Xi{x,z,z')e''^'^^-'^'''\{x-z')^{^) 



'\x\<^RJT. JT, 




IT, JT. JD{ri,z,z') 

X (x — 7) (x — ■^'Y Q^{x — j)gri{x — 'j')dxdzdz' . 



We begin estimating 

|Xi|<CA:i+l^l 



S JS JD(ri,z,z') 



Xi{x, z, z')\x - 7|l^l|x - 7'|l^le-^^(l^-^ll'+l^-^l")dx(izd/. 
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Now, using the estimate (15) with p = \f3\, a = 5k and s = |a; — 7I or |rE — 7'!, and continuing 
the estimate of /i, we have for a constant, C, independent of z and z' , 



l/'l^ 



ii.i<«i^i«^ 



W\ 



S JT, JD{ri,z,z') 



Xi{x,z,z')e-f'^\'^-"'\'+\'-"''\'Uxdzdz' 



JS Jt, JD{ri,z,z') 

<Ck[ [ [ xi{x, z, /)e-f (l--^l'+l--^'l')e-f l^-^'l' dxdzdz' 



<Ck I I e 8 



ifefli,_^/|2 



e|2-z'|2 



e-f (I-7P+I-7'P) dxdzdz'. 



D{ri,z,z') 



Here we have used the identity 



1, 



I |2 I I /|2 ol -|2 I -^1 /|2 

|x-7| +|x-7| =2|x-7| +-I7-7I ) 

and the fact that I7 — 7'! > ^9\z — z'\ on the support of xi- For the inner integral we can use 
Cauchy-Schwarz, together with the fact that D C ^irj; z) and D C J^(?/; z'), 

1/2 

■f {|x-7|2 + |z-7'P)d^ < ( / **= 

lD{ri,z,z') \Jq{ti;z) 

By a change of local coordinates we can show that 

e-f l^-^l' dx < Cfc(i-'^)/2. 



L 



n(ri;z') J 



(45) 

From this it follows that 
(46) 



/ 



n(rj;z) 






To show (45) for each z, we introduce local coordinates in the tubular neighborhood ^{rj; z) 
around the ray 7 in the following way: choose (smoothly in (s, z)) a normalized orthogonal basis 
ei(s, z), ..., erf_i(s, z) in the plane {x : (x — x(s; z))-p(s; z) = 0} with the origin at x(s; z). Since 
■s and z lie in compact sets, there will be an r/ > such that in the tube 0(77; z) the mapping 
from X to {s,y) defined by 

X = x(s; z) + yiei(s, z) -\ h y^-i • ed-i{s, z) 

will be a diffeomorphism depending smoothly on z, hence 

9x 



C(r,;2) ^|s|<Lo ^|j/|<»7 



9(y,s) 



dyds < Ck'^^-'^^/'^, 



where Lq is chosen such that |x(Lo; z)| > 5i? for all z G S. Letting A = sup^ ^/g^ |z — z'| < cxd 
be the diameter of S, we continue to estimate the (z, z')-integral left in (46): 










< Ck'^-'^, 
which concludes the estimate of /i. 
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6.2. Estimate of I2. In order to estimate I2 we use a version of the non-stationary phase lemma 

(see t22j). 

Lemma 6.4 (Non-stationary phase lemma). Suppose that u{x;C) G Cq°{Q x Z), where 0. and 
Z are compact sets and ip{x; Q E C°^{0) for some open neighborhood O ofO,xZ. IfVxip never 
vanishes in O, then for any iT = 0, 1, . . ., 



u{x- Qe^'^^^'^'^'^dx 






\X\<K 



where Ck is a constant independent of C,. 
We now define 



l2{z,z):= I l2{x,z,z')dx 

J\x\<5R 

= fei+l/^l / X2{x, z, z')e*'^(^'^'^')5(x; z')'^{^) 

JD(ri,z,z') 

X (x - 7)^(x - -f')I^Qr^{x - 'y)Qn{x - -f')dx. 



In this case, non-stationary phase Lemma 6.4 can be applied to I2 with C, = {z,z') G S x S to 
give, 



<CKk 



i+m-K 



(C(0,r?)|z-z'|)2i^-|Al7^(^,,_,,) 



<CKk 



i+m-K 



|A|<i^ 

|A|<i^l^ ^1 Vai+A2=A-^^(^'^'^') 

Ai<2/3 



1 



E 



dt 



y 



{x - -ff{x - j'fx2g'ggvQr, 

{x-jf{x-Yf 



e-^'^^dx 



dx^ [x2g'ggr,Q'rj] 



JfcV'rfa 



where g'^ = Qri{x — 7'), and we used the fact that \Vx'4^{x, z, z')\ > C{6, t])\z — z'\ on the support 
of X2, shown in Lemma 6.3 The constant Ck is independent of z and z' . By the bound (38) 



and since ^^ is uniformly smooth and x, z, z' vary in a compact set, \d^^ \x2g' gQ-qQ'r\ \ can be 
bounded by a constant independent of x, z and z' . We estimate the other term as follows. 



a 



Ai 



(x-7)'^(a;-7')^ <C ^ (x - 7)'^-^"(x - 7')^"^^^ 

Aii+Ai2=Ai 
Aii,Ai2</3 

<C Y^ |a;-7|l^l-l^iil |x-7 



'|l/3|-|Ai2| 



Aii+Ai2=Ai 
Aii,Ai2</3 



18 HAILIANG LIU\ JAMES RALSTON^, OLOF RUNBORG^*, AND NICOLAY M. TANUSHEV 

Now, using the same argument as for estimating Ii, we have 



Jd(v.z.z') L 



Aii+Ai2=Ai ^ " ' ' 
Aii,Ai2</3 

-|/3| + |All|-|g| + |Ai2l 

< C(A2)A; 2 I e 

'D(r,,z,z') 

< (^^(l-d)/2-|/3| + |Ai|/2 






f{(x-7)2+(x-7')2)(^a 



and consequently. 

/2 



|A|<i^ ' ' Ai+A2=A 

Ai<2/3 



< Ci,fe(3-<i)/2 ^ 



1 



|A|<X 



(|z-z'|\/fc)2^-|^l 



On the support of X2 the difference |z — z'| can be arbitrary small, in which case this estimate is 
not useful. However, it is easy to check that the estimate is true also for K = {), and I2 is thus 
bounded by the minimum of the K = {) and K > Q estimates. Therefore, 



< ck^^-^y^ 



mm 



1 



\A<^ \z - z'\^k 



2A'-|A| 



< CA;(3-'^)/2 ^ min 

\\\<K 



< Cfc(3-<i)/2 ^ 



1, 



1 



z-z'\Vkj 



<C- 



k(3-d)/2 



\\\<K 1 + (\z - z'\Vk] 1+ [\z- z'\Vk] 



Finally, letting A = sup^ ^'es 1^ — -z'l < 00 be the diameter of S, we compute 

/ '- 

'^'<^l+(\z-z'\y/k) 



SxE 



12(2, z') 



dzdz' <Ck^ 



K 



dzdz' 



3-d 

<Ck^ 



A 



1 



< Ck 



2-d 



l + {TVk)^ 
00 cd—2 



T^-^dT 



1 + e 



K 



d^ 



< Ck 



2-d 



if we take K > d — 1. This shows the I2 estimate, which proves claim (39). 

7. Another Superposition 

Specializing to p{x) = {x — y) ■ v^ one can also take the superposition with respect to v. We 
will carry this out for d = 3. Starting with an inversion formula for the Radon transform: 
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fix) = -^A ( l^^du ( y^_ ^_^/(y)d^J 1 , 

and noting that /^2 d-^ I(x~v)-u=o f^y^'^^v ^^nds to zero as |x| — )• oo when / G Cc(M'^), it follows 
that 




In other words 

as a distribution. Hence, ignoring p and the lower order term 



/ f{y)dAy\ =27r f /^^dy. 

J{x-y)-u=0 ) JR'^ F ~ y\ 

/ S{x ■ v)di 



2Tri e~^\^~y\'^ 

g{u,y,k)dv = -—- =def h{x;y,k), 

52 k \x-y\ 

and Jg2UGB{x;i'-,y)dv is a approximation to the outgoing solution to L„u = h satisfying the 
estimate (pil). 



Acknowledgments 

This article arose from work at the SQuaRE project "Gaussian beam superposition methods 
for high frequency wave propagation" supported by the American Institute of Mathematics 
(AIM), the authors acknowledge the support of AIM and the NSF. 

References 

[I] G. Ariel, B. Engquist, N. M. Tanushev, and R. Tsai. Gaussian beam decomposition of high frequency wave 
fields using expectation-maximization. J. Comput. Phys., 230(6):2303-2321, 2011. 

[2] V. M. Babic and V. S. Buldyrev. Short-Wavelength Diffraction Theory: Asymptotic Methods, volume 4 of 

Springer Series on Wave Phenomena. Springer- Verlag, 1991. 
[3] V. M. Babic and T. F. Pankratova. On discontinuities of Green's function of the wave equation with variable 

coefficient. Problemy Matem. Fizikt, 6, 1973. Leningrad University, Saint-Petersburg. 
[4] V. M. Babic and M. M. Popov. Gaussian summation method (review). Izv. Vyssh. Uchebn. Zaved. Radiofiz., 

32(12):1447-1466, 1989. 
[5] J.-D. Benamou, F. CoUino, and O. Runborg. Numerical microlocal analysis of harmonic wavefields. J. Comput. 

Phys., 199(2):717-741, 2004. 
[6] N. Bleistein. Mathematical methods for wave phenomena. Academic Press, INC. 1984. 
[7] S. Bougacha, J.-L. Akian, and R. Alexandre. Gaussian beams summation for the wave equation in a convex 

domain. Commun. Math. Set., 7(4):973-1008, 2009. 
[8] F. Castella and T. Jecko. Besov estimates in the high-frequency Helmholtz equation, for a non-trapping and 

C^ potential. J. Diff. Eq., 228(2) :440-485, 2006. 
[9] F. Castella, B. Perthame, and O. Runborg. High frequency limit of the Helmholtz equation II; Source on a 

general smooth manifold. Commun. Part. Diff. Eq., 27:607-651, 2002. 
[10] V. C. Cerveny, M. M. Popov, and I. Psencik. Computation of wave fields in inhomogeneous media — Gaussian 

beam approach. Ceophys. J. R. Astr. Soc, 70:109-128, 1982. 

[II] B. Engquist and O. Runborg. Computational high frequency wave propagation. Acta Numerica, 12:181-266, 
2003. 

[12] Y. A. Erlangga. Advances in iterative methods and precondit loners for the Helmholtz equation. Arch. Comput. 

Methods Eng., 15:37-66, 2008. 
[13] E. Faou and C. Lubich. A Poisson integrator for gaussian wavepacket dynamics. Computing and Visualization 

in Science, 9(2):45-55, 2006. 
[14] G. A. Hagedorn. Semiclassical quantum mechanics. I. The fi — >■ limit for coherent states. Comm. Math. 

Phys., 71(l):77-93, 1980. 
[15] E. ,1. Heller. Time-dependent approach to semiclassical dynamics. J. Chem. Phys., 62(4):1544-1555, 1975. 



20 HAILIANG LIU\ JAMES RALSTON^, OLOF RUNBORG^, AND NICOLAY M. TANUSHEV* 

[16] E. J. Heller. Frozen Gaussians: a very simple semiclassical approximation. J. Chem. Phys., 76(6):2923-2931, 

1981. 
[17] M. F. Herman and E. Kluk. A semiclassical justification for the use of non-spreading wavepackets in dynamics 

calculations. Chem. Phys., 91(l):27-34, 1984. 
[18] N. R. Hill. Gaussian beam migration. Geophysics, 55(11):1416-1428, 1990. 
[19] N. R. Hill. Prestack Gaussian beam depth migration. Geophysics, 66(4):1240-1250, 2001. 
[20] L. Hormander. Fourier integral operators. I. Acta Math., 127(l-2):79-183, 1971. 
[21] L. Hormander. On the existence and the regularity of solutions of linear pseudo-differential equations. 

L'Enseignement Mathematique, XVII:99-163, 1971. 
[22] L. Hormander, The Analysis of Linear Partial Differential Operators I: Distribution Theory and Fourier 

Analysis, Springer- Verlag, Berlin Heidelberg New York, 1983. 
[23] S. Jin, P. Markowich, and C. Sparber. Mathematical and computational models for semiclassical Schrodinger 

equations. Acta Numerica, pages 1-89, 2012. 
[24] S. Jin, H. Wu, and X. Yang. Gaussian beam methods for the Schrodinger equation in the semi-classical 

regime: Lagrangian and Eulerian formulations. Gommun. Math. Sci., 6:995-1020, 2008. 
[25] S. Jin, H. Wu, X. Yang, and Z. Y. Huang. Bloch decomposition-based Gaussian beam method for the 

Schrodinger equation with periodic potentials. J. Gomput. Phys., 229(13):4869-4883, 2010. 
[26] S. Jin, H. Wu and X. Yang, A Numerical Study of the Gaussian Beam Methods for One- Dimensional 

Schrodinger-Poisson Equations, J. Comp. Math., to appear. 
[27] A. P. Katchalov and M. M. Popov. Application of the method of summation of Gaussian beams for calculation 

of high-frequency wave fields. Sou. Phys. Dokl., 26:604-606, 1981. 
[28] A. Katchalov, Y. Kurylev and M. Lassas. Inverse boundary spectral problems. Chapman and Hall (2001) 
[29] J. Keller. Geometrical theory of diffraction. J. Opt. Soc. Amer, 52, 1962. 
[30] L. Klimes. Expansion of a high-frequency time-harmonic wavefield given on an initial surface into Gaussian 

beams. Geophys. J. R. astr. Soc, 79:105-118, 1984. 
[31] S. Leung and J. Qian. Eulerian Gaussian beams for Schrodinger equations in the semi-classical regime. J. 

Gomput. Phys., 228:2951-2977, 2009. 
[32] S. Leung, J. Qian, and R. Burridge. Eulerian Gaussian beams for high frequency wave propagation. Geo- 
physics, 72:SM61-SM76, 2007. 
[33] J. Lu and X. Yang. Convergence of frozen Gaussian approximation for high frequency wave propa- gation. 

Gomm. Pure Appl. Math., 65:759-789, 2012. 
[34] H. Liu and J. Ralston. Recovery of high frequency wave fields for the acoustic wave equation. Multiscale 

Model. Sim., 8(2):428-444, 2009. 
[35] H. Liu and J. Ralston. Recovery of high frequency wave fields from phase space-based measurements. Mul- 
tiscale Model. Sim., 8(2):622-644, 2010. 
[36] H. Liu, O. Runborg, and N. M. Tanushev. Error estimates for Gaussian beam superpositions. Math. Gomp., 

82:919-952, 2013. 
[37] A. Majda, and J. Ralston. An analogue of Weyl's theorem for unbounded domains, II. Duke Math. .Journal, 

45:183-196, 1978. 
[38] M. Motamed and O. Runborg. Taylor expansion and discretization errors in Gaussian beam superposition. 

Wave Motion, 2010. 
[39] J. Ralston. Gaussian beams and the propagation of singularities. In Studies m partial differential equations, 

volume 23 of MAA Stud. Math., pages 206-248. Math. Assoc. America, Washington, DC, 1982. 
[40] B. Perthame and L. Vega. Morrey-Campanato estimates for Helmholtz equations. Journal of Functional 

Analysis, 164:340-355, 1999. 
[41] M. M. Popov. A new method of computation of wave fields using Gaussian beams. Wave Motion, 4:85-97, 

1982. 
[42] J. Qian and L. Ying. Fast Gaussian wavepacket transforms and Gaussian beams for the Schrodinger equation. 

J. Gomput. Phys., 229:7848-7873, 2010. 
[43] V. Rousse and T. Swart. A mathematical justification for the Herman-Kluk propagator. Gomm. Math. Phys., 

286(2):725-750, 2009. 
[44] O. Runborg. Mathematical models and numerical methods for high frequency waves. Gommun. Gomput. 

Phys., 2:827-880, 2007. 
[45] N. M. Tanushev. Superpositions and higher order Gaussian beams. Gommun. Math. Sci., 6(2):449-475, 2008. 
[46] N. M. Tanushev, J. Qian, and J. V. Ralston. Mountain waves and Gaussian beams. Multiscale Model. SimuL, 

6(2):688-709, 2007. 
[47] N. M. Tanushev, B. Engquist, and R. Tsai. Gaussian beam decomposition of high frequency wave fields. J. 

Gomput. Phys., 228(23):8856-8871, 2009. 



21 

[48] B. Vainberg. On short-wave asymptotic behaviour of solutions to steady-state problems and the asymptotic 
behaviour as i — > oo of solutions of time-dependent problems. Uspekhi (Russian Math. Surveys), 30(2): 1-58, 
1975. 

[49] B. R. Vainberg. Asymptotic Methods in Equations of Mathematics Physics, Gordon and Breach (1989) 

Appendix A. Form of the Green's Function 

Let G\{x) be the free space Green's function for the Helmholtz equation at complex valued 
wave number A = |A|/3 where (5 is complex number with |/3| = 1 and 9/3 > 0. The Green's 
function has the following properties, 

(47) Ga(x) = 0(e-^'=l^l|x|^), drGx{x) - iXGx{x) = 0{\x\'-^), r = |x| ^ oo. 

The dependence on \k\ can be scaled out and by rotational invariance we can write Gx{x) = 
\X\'^-^Gp{\Xx\) where G^{x) = Gfi{\x\). Then, if 



[pr) 2 
the complex valued function Wj^ will satsify the following ODE for r > 0, 

2 



(48) 



Cd 



wl{r) + 2ipw'p{r)-^wp{r)=0, 



Cd 



d-2 



This follows from applying the Helmholtz operator in d dimensions to Gp away from x = (with 

r = \x\), 

= ^Gp{x) + P^Gp{x) = ^G^ir) + ^i-G'^(r) + /3^Gp{r) 



oipr 



(/3r) 



■JZT (w"fi{r) + 2ifiw'p{r) 



r dr 
{d-l){d-2>)wp{r] 



After differentiating (|48|) p times we get 



(49) 



w 



(P+2), 



-,(P+1)/ 



3=0 



jwf{r)r 



-2-p+j 



0, 



for some coefficients dpj. Prom the left property in (47) it follows that |it;^(r)| < Bq for some 
bound Bq and r > 1. Moreover, the right property (the radiation condition) implies that 



a —r {d — l)wp/2r as r — )■ oo. It then follows by induction on (49) that wF {r) — )■ for all 



w 
p> 1. 

We now claim that there are bounds Bp, independent of r, such that \rPiI)o' {r)\ < Bp for 
r > 1. We just saw that this is true for p = and we make the induction hypothesis that it is 
true for j = 0, . . . ,p. Then from (49), 

^-2rS/3 
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when r > 1, where i^p+i = X]f=o l^pj^jl- Since Wg (r) — ;• as r — > oo and 9/3 > 



-(p+i)/ \ 
Wg (r) 



„2r3/3 



ds 



,'^ifis^(P+l)(s)ds 
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where -Bp+i = B'p^^/{p + 1). This shows the claim. 
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We conclude that 

i\\x\ 

Ga(x) = \\\''-^Gp{\\x\) = ^^^w{x-\), w{x-\) = \\\'^p'-^wp{\\x\), 

\x\ 2 



and for any multi-index a, 
|a>(x;A)|<C|A|^j; 

<C((5)|A|^, 
when |a;| > 6 and |A| > 1/5. 



dr'J 



w,3{\\\r) 



W\ 



W\ 



r=\x\ 



lAi'^EhH'H^N) =i^i^E^^- 



j=0 



i=o 



